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Abstract 

The spectral problem associated with the linearization about solitary waves of spinor systems 
or optical coupled mode equations supporting gap solitons is formulated in terms of the Evans 
function, a complex analytic function whose zeros correspond to eigenvalues. These problems 
may exhibit oscillatory instabilities where eigenvalues detach from the edges of the continuous 
spectrum, so called edge bifurcations. A numerical framework, based on a fast robust shooting 
algorithm using exterior algebra is described. The complete algorithm is robust in the sense 
that it does not produce spurious unstable eigenvalues. The algorithm allows to locate exactly 
where the unstable discrete eigenvalues detach from the continuous spectrum. Moreover, the 
algorithm allows for stable shooting along multi-dimensional stable and unstable manifolds. 
The method is illustrated by computing the stability and instability of gap solitary waves of 
a coupled mode model. 
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1 Introduction 



Transmitting information efficiently across long optical waveguides is a big challenge in telecom- 
munications. Gap solitons are potential candidates to achieve this goal. A gap in the linear 
spectrum allows solitons in spinor-like systems to propagate without losing energy due to a res- 
onant interaction with linear waves jlll I18j . For example, an optical fiber with a periodically 
varying refractive index supports gap solitons. The gap here is created by Bragg reflection and 
resonance of waves along the grating. 

Before its application to optical waveguides and transmission of optical pulses, gap solitons 
have been studied in the context of spinor field equations in elementary particle physics ^1 120] 
and in condensed matter physics 0. There is an extensive literature on the existence of gap 
solitons, but the issue of their stability, which is of paramount practical importance, is still an 
open question for many systems. 

Gap solitons were long believed to be stable. This was conjectured on the grounds of computer 
simulations in (restricted) parameter regimes. Only recently Barashenkov, Pelinovsky 
& Zemlyanaya showed analytically by using a perturbation theory, and Barashenkov & 
Zemlyanaya verified numerically that, in a particular optical system, gap solitons can un- 
dergo oscillatory instabilities where eigenvalues detach from the edges of the continuous spectrum 
(edge bifurcations). In Kapitula & Sandstede Evans functions are used to detect analyt- 
ically the onset of an oscillatory instability (edge bifurcation) in near integrable systems. From an 
analytical point of view a difficulty in proving stability/instability is that the energy in these sys- 
tems is neither bounded from above nor from below. Usual techniques such as energy-momentum 
methods or methods involving Lyapunov functionals are therefore bound to fail. 

An often used numerical approach to eigenvalue problems is to discretize the spectral problem 
on the truncated domain x £ [—Loo, Lao] (with Lqo ^ 1) using finite differences, collocation or 
spectral methods, reducing it to a very large matrix eigenvalue problem. There are two central 
difficulties with this approach. Firstly, in general the exact asymptotic boundary conditions 
at X = ±Loo depend on the eigenvalue A in a nonlinear way, and so application of the exact 
asymptotic boundary conditions changes the problem to a matrix eigenvalue problem, which is 
nonlinear in the parameter. So matrix eigenvalue solvers can no longer be used. In the above 
papers artificial boundary conditions such as Dirichlet or periodic boundary conditions, were 
applied, in order to retain linearity in the spectral parameter. For a detailed discussion on artificial 
boundary conditions, see Sandstede & Scheel Secondly, the approximate boundary 

conditions lead often to spurious discrete eigenvalues generated from the fractured continuous 
spectrum. If the continuous spectrum is strongly stable (that is, the continuous spectrum is 
stable and there is a gap between the continuous spectrum and the imaginary axis) this does 
not normally generate spurious unstable eigenvalues. However, if the continuous spectrum lies 
on the imaginary axis (as often happens with gap solitons due to the Hamiltonian character 
of the underlying system), spurious eigenvalues may be emitted into the unstable half plane. 
Indeed, Barashenkov & Zemlyanaya give an extreme example, where a large number of 
spurious unstable eigenvalues are generated by the matrix discretization (see Figure 1 in [5]). 
These problems prohibit their method to be used as a robust scheme to check stability when no 
theoretical results are given to guide the experiments. 

A robust and stable numerical scheme is still an open problem. In this paper we present a nu- 
merical method based on the Evans function and exterior algebra which does not exhibit spurious 
unstable eigenvalues. The Evans function is a complex analytic function whose zeros correspond 
to eigenvalues of the spectral problem associated with the linearization about a solitary wave so- 
lution. The Evans function was first introduced by EvANS and generalized by Alexander, 
Gardner & Jones ,2 . To define the Evans function, one writes the wave equations as a system 
of first order real equations with respect to the spatial variable x , such that one gets a system 



of the form Z^. = -F(Z, Z^, Z^^, . . .) . The study of hnear stabihty of a sohtary wave solution Z(x) 
involves a linearisation about Z by writing the basic solutions as Z(x) = Z(x) + u{x) e^^ . This 
will lead to a linearized problem of the form 



where A E C is the spectral parameter and A(x, A) is a matrix in C"^" , whose limit for x ±00 
exists. The solitary wave solution Z of a partial differential equation is linearly unstable if for 
a spectral parameter A with 3f?(A) > there exists an associated perturbation u(x) , which is 
bounded for all x . An oscillatory instability (edge bifurcation) does happen when the continuous 
spectrum is on the imaginary axis and one of the bounded eigenfunctions of the continuous 
spectrum develops into an exponentially decaying eigenfunction with 3f?(A) > 0. 

To verify if such bounded eigenfunction u{x) exists for a given value of A with 5R(A) > 0, 
it is checked if the unstable manifold at x = — 00 and the stable manifolds at x = 00 have a 
non-trivial intersection. If this is the case, the solitary wave solution Z is unstable with positive 
growth rate 3f?(A) > 0. To check the transversality condition the Evans function is used. Hence 
the Evans function can be viewed as a Melnikov function or a Wronskian determinant. 

Crucial to the initial construction of the Evans function is the distribution of the eigenvalues 
of the 'system at infinity', that is, the matrix A±oo(A) which is associated with the limit as 
X — > ±00 of A(x, A) . It is assumed that no eigenvalues are on the imaginary axis ^ and that the 
number of negative eigenvalues is constant for A € A , where A is a simply-connected subset of 
C. Let k be the number of negative eigenvalues of A-i-oo(A) for A G A. Note that in the case of 
different asymptotic behaviour at x = ±00 as for example in the case of fronts, k can be different 
at X = ±00 in general. 

This suggests a naive approach by which one may follow the stable/unstable manifolds at 
X = ±00 with a standard shooting method and check their intersection using the Evans function. 
This may be done indeed if the dimension of these manifolds is 1. Otherwise, any integration 
scheme will inevitably just be attracted by the eigendirection corresponding to the most unstable 
eigenvalue. However, for most systems, the dimensional of the stable/unstable manifold will be 
larger than 1. In this paper, we will consider a model equation with stable/unstable manifolds of 
dimension 2. To keep the eigendirections orthogonal in the course of the numerical integration, 
one may employ a Gram-Schmidt orthogonalization method. However, this is a non-analytic 
procedure, which will eventually backfire in the case of oscillatory instabilities where we expect 
zero's of the Evans function in the complex plane. Indeed, to locate those complex eigenvalues, 
Cauchy's Theorem (argument principle) will be employed and thus an analytical method is crucial. 
In our numerical method, we will use exterior algebra, which allows for an analytical calculation 
of the Evans function. 

The numerical method builds on the work in |H1 llOj , where a numerical algorithm is given for 
the calculation of the Evans function and applications are give to the stability and instability in 
a fifth-order KdV equation. In this paper, we will show that a similar algorithm can be used to 
determine oscillatory instabilities (edge bifurcations). The algorithm does not exhibit spurious 
eigenvalues and the exact asymptotic boundary conditions are built into the definition of the Evans 
function in an analytic way. The analyticity can then be utilized to apply Cauchy's principle value 
theorem to study stability/instability. An important feature of the numerical method is that it 
involves the use of exterior algebra to describe the system on a higher dimensional space in which 
a simple shooting method can be employed. A similar idea is used in Brin 7^ and Brin &; 
ZuMBRUN j8j for dealing with instabilities in viscous fluid flows. 

^The assumption on the hyperbolicity of A±oo(A) can be weakened (see |13II14| '). 



Ua, = A(x,A)u, U G C 




To illustrate our method we will consider the following coupled mode model 

= i{ut + u,j;) + v + + /9|np)n 
= i{vt -Vx)+u + (liip + p\v\'^)v 

This is a model to describe optical pulses in waveguides which have been grated so the refractive 
index is varying periodically. In the case p = this equation is known in field-theory as the 
massive Thirring model and was shown to be completely integrable (see for example ^ ) . In a 
nonlinear optics context, one has p = 1/2 in periodic Kerr media jJJ, but in other media p may 
range from up to infinity |^. The equation H1.2j) has also been studied by Barashenkov, 
Pelinovsky & Zemlyanaya Barashenkov & Zemlyanaya 5^ and in a shghtly modified 
form by Kapitula & Sandstede ^Sj. In |^ |S] a heuristic perturbation analysis is used to 
analyse the onset of oscillatory instabilities and numerical study is used to give a more complete 
picture. The numerical study encountered serious problems as discussed above. In |15| . a relation 
between the Evans function and the inverse scattering formalism is established for integrable 
systems. This forms the basis of a rigorous perturbation analysis for perturbations of the massive 
Thirring model. 

The coupled mode model (|1.2)) is chosen for illustration purposes, since the work in the 
previous papers allows us to illustrate the advantages of our method. We stress though that 
the method we present is general and can be applied to gap solitons in other systems as well. 
Moreover, there is no need for the system to be related to an integrable system. 



(1.2) 



2 The numerical Evans function for the coupled mode mode 

The solutions and dynamics of (|1.2)) are best described by splitting off the real and imaginary 
part of the fields n, v . We shall write for a solution u = Qi + iPi and v = Q2 + iP2 and collect 
the information in a single real solution vector Z = {Qi, Q2, Pi, P2) ■ 

2.1 The model equations and its solutions 

We introduce the Lorentz transformation X = [x — Vt)/\/\ — and T = {t — Vx)/\/l — . 
In these boosted variables the model system (|1.2|) may be written in a (quasi-)multisymplectic 
framework as 

E[MZt + ICZx] = VS{Z) (2.1) 
where Z = (Qi, Q2, A, ^2) , 

S{Z) = ^ [Q1Q2 + P1P2 + (Ql + Pfml + P|)] + f [iQl + P?f + (Ql + P2f] , 

I e-y 

E = I ^ I , El = I I , M = I " I and 

with V = tanh y and the Pauli matrices ctq 3 are defined as 



Note that if p = , the system is invariant under this Lorentz transformation. 








The reason for introducing this formahsm is that a multisymplectic formulation allows for a 
systematic linearization and also sheds light on conservation properties within the system. We 
note that (|2.1|) is equivariant under action of the continuous symmetry group S0{2) acting on 

represented by 



cos{ip) ao — sin(-(/;)(To 
sin(^/;)fTo cos(^/;) (Jq 



(2.2) 



since [G^,M] = [G^,K] = [G^,E] = and S{G^Z) = S{Z) for any ip . According to Noether's 
Theorem, there is a conservation law associated with this continuous symmetry, namely Pt + 
Qx = with P and Q determined by 



VP(Z) 



_d_ 



G^(Z) and VQ(Z) 



i>=o 



d_ 



(2.3) 



hence 



P(Z) = Y,Ql + Pf and Q{Z) = Y,{-^r^\Ql + Pf 



1=1 



i=l 



Note that in the original system (|1.2|) . this symmetry shows up as an equivariance of the system 
under simultaneous phaseshifts of u and v, i.e., u i— > ue^^ and v ^ ve 



Going to a frame moving with the symmetry group and writing the solutions as Z{X,T) 
G(p(x)-otZ(X, T) , the equations (|2.1|) become 



E[MZt + KZx - nVP{Z) + ^x^QiZ)] = V5(Z) , 



(2.4) 



where we dropped the tildes. Time independent solutions in the moving frame were found by 
AcEVES & Wabnitz P to be 



Z{X) = aE'2 




where 



= cos (0 < 6* < vr) 
1 

a 



( Wr{X)\ 
Wr{X) 

Wi{X) 
\Wi{X)j 



ip{X) = 2a^psinh(2y) arctan {tanh[(sin(6l))X] tan |} 

sm{e) 



(2.5) 



and W{X) = Wr{X) + iWi{X) 



cosh ((sin 6i)X- 16*72) ■ 



"v/l + pcosh 2y 

Note that there cannot be any gap solitons for \Q\ > 1 as then there is no gap in the linear 
spectrum. At the upper edge of the gap, at 6^0 (0^1), the solutions approaches the small- 
amplitude nonlinear Schrodinger soliton VF(Ar) = 6sech{9{X — [/2)) . At the lower edge, at ^ ^ tt 
(0 — > —1), the gap soliton has a finite amplitude and decays algebraically, M^(X) = i/(X + i/2). 
These two limits are referred to as "low intensity" and "high intensity" limits, respectively 



2.2 The linearized problem 

To study the stability of the solution (|2.5() , we linearize around this solution in the moving frame 
and use a spectral Ansatz 

Z(X,r) = (Z{X) + u{X)e^^^ . 

We obtain the dynamical system 

ux = A{X,X)u, ugC^ (2.6) 

with 

A(X, A) = [E-'^D^S{Z) + QD^P{Z) - ipxD^Q{Z) - AM] (2.7) 
where denotes the Hessian. 



2.3 Asymptotic properties of the linearised system 

As described in the Introduction, the eigenvalues and eigenfunctions of the asymptotic matrix 
Aoo(A) 

A(x, A) has the asymptotic property that 



lim A(x, A) are crucial for the numerical Evans function approach. The matrix 

X— >itoo 
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-A 







-ey^ 
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e~y 
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Aoo(A) = 


lim A(x, A) = 












X— >iboo 
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ey 


-A 









v- 


-e-y 


-n 





^ J 


The characteristic polynomial of 


A 


3o(A) 


is 






A(^,A) = 


det[^I- Aoo(A)] 






-A2- 


1 + 


f + 



(2. 



- 4:^ 2^ A^ 

Thus the eigenvalues /i of the asymptotic matrix satisfy 

- A^ - 1 + = ±2il^A or ^ 1 + (A ± i^f . 



(2.9) 



(2.10) 



The continuous spectrum is found by setting TZ{fi) = or /i = i/c . A short calculation gives 
that there are four branches of continuous spectrum on the imaginary A-axis: 



X±{k) = + k'^± \n\) and -A±(/c), k£R. 

The end points of the continuous spectrum are at ±i(l — and ±i(l + ■ For all values 
of A not in the continuous spectrum, there are 2 eigenvalues n with positive real part and 2 
eigenvalues fi with negative real part. Thus outside the continuous spectrum, there is a 2 : 2 
splitting, that is the dimension of the stable/unstable manifolds is 2 at both limits of x = ±00. 

With ()2.1U|) . it is easy to determine explicit expressions for the eigenvalues and eigenvectors. 
The eigenvalues are 



/i±(A) = ±v(A - i|oi)2 + 1 = ±^A - mi + 1) - m\ - 1), 



^±(A) = ±V(A + W + 1 = ±\/A + m\ + 1) + Kl^^l - !)■ 



where the square roots are defined as follows: 

^= y|i|ei'^''g(^)/2, where - vr < arg(2) < vr and z = X ±i{\n\ ±1). 

Hence the cuts in the complex plane, associated with this definition of the square roots, start at 
the end points of the continuous spectrum and continue in the left half of the complex A-plane. 
Note that this definition implies that all square roots have positive real parts. 

In this model system it is easy to find explicit expressions for the eigenvectors. The eigenvector 
with the eigenvalue /i^ for the matrix Aoo(A) is given by 

vec±(A) = (F±,sgn(J7)ie-^-sgn(J])iF±,e-?^)^, where F± = A - ^± - i|f]|, 
Yec^{X) = {F^,-sgn{n)ie-y,sgn{n)iF^,e-y f, where = X - fi^ + 

The eigenvector with the eigenvalue for the adjoint matrix (Aoo(A))^ is given by 

advec^(A) = (Fm,sgn(0)ie^,-sgn(0)iFm,e^)'^, where F^ = X - fj.^ - 
advec^{X) = (F^,-sgn{Q)iey,sgn{n)(F^,ey)'^, where F^ = X - fi^ + i\n\. 

One might also use a numerical routine to solve the eigenvalue problem of Aoo(A) . However, 
one has to be very careful near the end points of the continuous spectrum. At these points two 
eigenvalues collide, hence it will be tricky to find numerically the 'correct' eigenvalue. If we fail 
to do so, this will result in a change of orientation in the Evans function 1)2.14(1 . This causes 
unwanted zero's of the Evans function and is reminiscent of the unstable spurious eigenvalues 
observed when the system (|'2.6() is solved by some form of discretization as done by Barashenkov 
& Zemlyanaya i5j. One can construct a numerical routine to follow the correct eigenfunction, 
but if the analytical expressions are readily available (as they are here) it is much easier. 



2.4 The Evans function and the formulation in the exterior algebra 

For 3?(A) > 0, the system ()2.()|) and the properties of the system at infinity, Aoo(A) , are in stan- 
dard form for the dynamical systems formulation of the spectral problem proposed by EvANS |12j 
and generalized by Alexander, Gardner & Jones ;2,. A value of A G A is an eigenvalue if 
the 2 -dimensional space of solutions which decays as x — > — oo and the the (4 — 2) -dimensional 
space of solutions which do not grow exponentially as x — > +oo have a nontrivial intersection. 
The Evans function is an analytic function which gives a zero if such a nontrivial intersection 
exists. To obtain an analytic description of the 2 -dimensional space of solutions of ()2.6j) which 
do not grow exponentially as x ^ +oo, we will use that the system (|2.6() induces a dynamical 
system on the wedge-space A^(C^) . This is a space of dimension (2) = 6 . To define the Evans 
function, the induced dynamics on this wedge-space A (C^) will be used. 
The induced system can be written as 

U, = A(2)(x)U, UeA'(C^)- (2.11) 

Here the linear operator (matrix) A^^^ is defined on a decomposable 2-form ui A U2 , Uj G C^, 
as 

A(2)(ui a U2) := (Aui) A U2 + ui A (Aug) (2.12) 

and it extends by linearity to the non-decomposable elements in /\^{C^) . This construction can 
be carried out in a coordinate free way and can be generalised to {k):{n — k) splittings in n- 
dimensional dynamical systems. General aspects of the numerical implementation of this theory 
can be found in Allen & Bridges 3_. 



Since the induced matrix A(^)(x,A) inherits the differentiabihty and analyticity of A(x,A), 
the hmiting matrices wih exist, 



Ag\X)= lim A^'\x,X). 

(2) 

The set of eigenvalues of the matrix Aoo (A) consists of aU possible sums of 2 eigenvalues of Aoo(A) 
(this is an exercise in multi-linear algebra, see Marcus H^)- Therefore, for K(A) > 0, there is an 
eigenvalue of A^xd (A) , denoted by o"_|_(A), which is the sum of the 2 eigenvalues of Aoo (A) with 
negative real part, i.e., o"+(A) = — (/x,^(A) + iJ,^{\)) (note that the subscript " + " in cr+(A) refers 
to exponentially decaying behaviour at +oo). Moreover this eigenvalue is simple, an analytic 
function of A and has real part strictly less than any other eigenvalue of A^'*(A) . Similarly, there 
is an eigenvalue cr_(A), which is the sum of the 2 eigenvalues of Aoo(A) with non-negative real 
part, i.e., o"_(A) = fJ-^W + f^pWi ^-('^) is simple, an analytic function of A, and has real 

(2) 

part strictly greater than any other eigenvalue of AJ^ (A). Note that cr_(A) = — (j+(A) in this 
example. 

Let C^(-^) be the eigenvectors associated with a±{X) , defined by 

A(^)(A)C+(A) = a+(A)C+(A) and A(^)(A)C-(A) = a^{X)CW ■ (2-13) 

These vectors can always be constructed in an analytic way (see jlTH). From section ESI it follows 
that 

C±(A) = vec±(A)Avec±(A). 
This implies 

C± = (F± + F±, -2e^F±F±, sgn(17)i(F± - sgn(17)i(F± - F±), -2e-^ F± + , 

The solution U^(x,A) is the solution of the linearised system (|2.1H) with the property that 
lim^j^-l-oo e~'^*^'^''^U^(x, A) = ^^(A). In this example, dimension of the unstable manifolds at 
X = oo and x = — oo are the same. For the construction of the general case see p| llOj. 

Note that the solutions U^(2;,A) are analytic expressions which represent the space of solu- 
tions which decay as x ^ ±00 in the original system (|2.6)) . With this the Evans function can be 
defined as 

^(A) =e-^o'^M)'^^U-(x,A) AU+(x,A), AG A, (2.14) 

where A is the wedge product and 

r(x,A) =Tr(A(x,A)). (2.15) 

For the case of the coupled mode equation <\l.'2\i this expression simplifies, since Tr(A(x, A)) = 
0, see (|2.7|) . The Evans function, as defined above, can be extended across the continuous 
spectrum with some cuts in the complex plane, see Gardner & Zumbrun |13j and Kapitula 
& Sandstede 14.. 

Next we will give an equivalent description of the Evans function, using the adjoint system. 
The adjoint system of (|2.11|) is 

-T 



V, = -[A(2)(x)] V. (2.16) 

The dimension of the unstable manifold of this adjoint system is equal to the dimension of the 
stable manifold of the linearised system (|2.11|) and its most unstable eigenvalue is —0"+ . Using 



the Hodge-star operator, we can relate the most unstable solution at x = — oo of the adjoint 
system with the most unstable solution U~ of the linearised system at x = — oo. Details can 
be found in Ell- ^o formulate the alternative description of the Evans function, we have to 
construct an inner product on /\^(C^) . Let 

U = ui A U2 and V = vi A V2 , Uj, Vj G C ^ , V z, j = 1, 2 , 

be any decomposable 2 -forms. The inner product of U and V is defined by 

(ui,Vi) (U1,V2) 
(U2,vi) (U2,V2) 

where (•,-)4 is the complex inner product in C^. Since every element in /\^(C^) is a sum 
of decomposable elements, this definition extends by linearity to any 2-form in /\^(C^) . An 
equivalent definition of the Evans function (|2.14jl is given by the following readily computable 
expression 

i^(A) = [V-(0,A),U+(0,A)l2, (2.17) 

where V~(x,A) is the most unstable solution at x = —oo of the adjoint system 1)2. 16(1 . From 
__ J' 

Section 12.31 it follows that the eigenvector of the adjoint matrix — (A^^)) for the eigenvalue 
-a> = Hp + is 

V' = [f+ + F+, -2e?^Fp+F+, sgn(0)i(F+ - F/), sgn(J7)i(F+ - F+), -2e-^ F+ + . 

Hence the solution V~(x,A) is the solution of the adjoint linearised system ()2.16|) with the 
asymptotic behaviour lim e'^+('*')^V~(x, A) = r/~(A). The generalisation of this definition for 

general splittings and more details can be found in |6l ITfl]. 

For the numerical implementation, we will need a basis for /\^(C'^) , and the above construc- 
tion assures that any basis will do. Therefore there is no loss of generality in assuming that the 
bases chosen are the standard ones. Starting with the standard basis for C^, and volume form 
V = ei A • • • A 64 , let ai, . . . , ag be the induced orthonormal basis on /\^(C^) . Using a standard 
lexical ordering, this basis can be taken to be 



[U,Vl2 :=det 



ai = eiAe2, a2 = eiAe3, a3 = eiAe4, 

(2.18) 

a4 = e2Ae3, as = e2 A 64, a6 = e3Ae4. 
Any U G A (C^) 

can be expressed as U = Yl^=i ^j^j ■ Since the basis elements a^ are orthogo- 
nal and the inner product |-, •]2 on /\^(C^) is equivalent to the inner product (•, •)e on C^, the 
expression (j2.17|) for the Evans function can be expressed in the equivalent form 

F(A) = (V-(0,A),U+(0,A))6. (2.19) 

It is this form of the Evans function we implement in our numerical algorithm. 
The matrix A^^) : A^(C^) ^ A^(C^) 

can be associated with a complex 6x6 matrix with 

entries 

{A(2)},^^. = [a„A(2)aj]2, i,j = l,...,6, (2.20) 



where, for any decomposable U = ui A U2 € /\^(C^) , A^^^U := Aui A U2 + ui A Au2. Let A 
be an arbitrary 4x4 matrix with complex entries aij,i,j = 1, • • • ,4, then, with respect to the 
basis (|2.18jl . A^'^) takes the explicit form 



aii + a22 


023 


024 


-Oi3 


-Oi4 





032 


011 + 033 


034 


012 





-Ol4 


042 


043 


flu +044 





012 


Ol3 




021 





022 + 033 


034 


-024 


-041 





021 


043 


022 + 044 


023 





-041 


031 


-042 


032 


O33+O44 



Details for these constructions in more general systems can be found in |^. 



2.5 Integration scheme 

The subtle nature of the oscillatory instability requires a high order integration scheme for the 
numerical integration of the linearised and the adjoint system. Instead of using the second order 
Gauss-Legendre Runge-Kutta method, i.e. the implicit midpoint rule, as in Bridges, Derks & 
GOTTWALDfTDI, we employ here a fourth order Gauss-Legendre scheme. We solve in 

Un+i = u" + ^Ax(Ki + K2), (2.21) 

where Ax is the spatial step size, sub- and superscripts n denote the spatial discretization and 
Ki_2 are implicitly defined by 

Ki = a(2)(x„ + (1 + ^)Ax) X (U'^ + iAxKi + (i + ^)AxK2) 
2 b 4 4 6 

K2 = A(2)(x„ + (i-^)Ax)x(U- + iAxK2 + (^-^)AxKi). (2.22) 
2 b 4 4 6 

In practice we solve (j2.22|) for and then subsequently we can solve (|2.21|) for U"^-'- . Since 

all equations are linear, the implicit form of H2.22() can be cast in an explicit form. 

The procedure for the numerical calculations is as follows. As explained in Section [2. 41 it is 
sufficient to restrict the shooting algorithm to /\^(C'^) . As a starting vector for the shooting 

algorithm we use the eigenvectors of Aoo (A) in the far-field (see Sections 12. ,31 and 12. 4|) . For the 
integration of the linearised system starting at x = +L00 (with L^o ^ 1) the starting vectors 
for each A are the eigenvectors C'''(A) related to the eigenvalues with the largest negative real 
part; for the integration of the adjoint system starting at x = — Lqo the starting vectors are 
the eigenvectors t]~{X) related to the eigenvalues with the largest positive real part. We have 
build in an analytical normalisation process such that the eigenvectors are normalized so that 
(r/~(A), C"''(A))6 = 1, for large values of A. (Note that this normalisation is not used for values of 
A < 2.) 

To calculate the Evans function, the linearised equation on /\^{C^) 

= [A(2)(x, A) - c7+(A)Irf]U+ , U+(x, A)|^^^^ = C+(A) , (2.23) 

is integrated from x = L^o to x = , where the scaling 

U+(x, A) = e-"+W^ U+(x, A) (2.24) 



ensures that any numerical errors due to the exponential growth are removed and U^(j;, A) = 
U+(x,A)|^_Q is bounded. An alternative to this scaling is to impose a renormalization of the 

vectors during or at the end of the integration, with for example |U"''(0,A)| = 1, but such a 
scaling does not preserve analyticity. 
For X <0, the adjoint equation 

^V- = [-A(2)(x,A)^ + MA)Id]V- , V-(x,A)|^_^^=r?-(A), (2.25) 

is integrated from x = —L^ to x = 0, also using the implicit midpoint rule, where again we 
introduce a rescaling 

V- {x, A) = e^+(^^' V- (x. A) (2.26) 

to remove the exponential growth. 

At X = 0, the computed Evans function is 

E{\) = (V-(0, A), U+(0, A))6 = (V-(0, A), U+(0, A))6 . (2.27) 



3 Numerical results and discussion 

In this Section we show results of our algorithm for the detection of oscillatory instabilities in the 
perturbed massive Thirring model H1.2() . We do not attempt here to present a thorough numerical 
analysis of the bifurcation scenarios of (|1.2() . The reader is referred to [5]. Instead our objective 
here is the presentation of a numerical algorithm which uses the Evans function as a numerical 
diagnostic tool for analysing edge bifurcations. Therefore we limit ourselves to illustrating several 
features of our numerical method. In Figure we have sketched a cartoon which indicates where 
in parameter space, i.e., in the 0-/j-plane, the numerical analysis takes place. This cartoon is a 
guidance to help to locate the upcoming plethora of figures and how these figures are related to 
bifurcations and instabilities. 
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Figure 1: Cartoon to illustrate where in parameter space the numerical analysis takes place. 

We study instability by computing the Evans function E[\) as defined in (|2.19|) . The analyt- 
icity of the Evans function for 5R(A) > allows one to detect oscillatory instabilities, i.e. complex 
roots of the Evans function, by means of Cauchy's theorem. The winding number of a closed 
curve in the A -plane tells us about the number of unstable eigenvalues. In all our calculations, 
we compute the complex Evans function E{\) , while varying the spectral parameter A = iA on 
the imaginary axis or varying the spectral parameter parallel to the imaginary axis with a (small) 
offset, explicitly A = off-|-iA, where "off' is the offset. The normalisation of the Evans function is 



chosen such that the Evans function converges to 1 for A large. Hence the closed curve is formed 
by connecting the endpoints of the imaginary axis via the half-circle with infinite radius. On this 
half-circle, the Evans function will always have the value 1. 

Since the system ()1.2() has translational and rotational symmetry, the Evans function will 
have a fourth order zero at A = 0. This means that in the vicinity of A = 0, the Evans function 
scales as E{X) ~ A^. Hence even an offset of only 10~^ yields that E{X) is of the order of 
10~^^ , making the calculations meaningless. In order to avoid this problem, we use an offset 
of at least 5 • 10~^ to analyse the Evans function near A = 0. Since we are interested here in 
oscillatory instabilities which occur at the edges of or within the continuous spectrum, and not in 
translational instabilities where eigenvalues emanate from A = 0, the offset near the A = does 
not affect the detection of instabilities. 




Figure 2: Real (continuous red line) and imaginary (dashed blue line) part of the Evans function 
E(X) as function of the spectral parameter A with off = 0, for V = 0.9 and 9 = OAn 

(a) The integrable case p = , The inset is a blow-up of the Evans function near the edge of the 
lower branch (X- ) of the continuous spectrum. 

(b) : p = 0.1. The inset is a blow-up of the Evans function and shows the discrete eigenvalue 
where E(X) = . The endpoints of the continuous spectrum occur at the cusps of the Evans 
function. 

In the integrable case {p = 0), the linearised massive Thirring model ()1.2|) has only neutral 
and continuous eigenvalues I17j and the solitary wave is stable for any < ^ < vr . The 
Evans function has zeros at the end points of the branches of continuous spectrum if p = . In 
Figure 12^, we used our algorithm to calculate the Evans function for p = , while the spectral 
parameter A is on the imaginary axis (off=0). This illustrates the zeros of the Evans function at 
the end points of the branches of continuous spectrum. As follows from Section 1231 the edges of 
the continuous spectrum are located at A-|- = i{\i^\ -\- 1) (upper branch) and at A_ = i{l — \^}\) 
(lower branch). At these endpoints, the real or the imaginary part of the Evans function exibits a 
cusp, as can be seen in Figure|2K. These cusps illustrate the non-analyticity of the Evans function 
at those points. Note that the cusps and the associated non-analyticity of the Evans function at 
the endpoints of the continuous spectrum occur for all p > and all < < vr , as can be seen 
from the following figures. 

For < ^ < 7r/2 and p > 0, the linearised perturbed massive Thirring model (|1.2j) has a 
discrete eigenvalue, which is located on the imaginary axis in the gap between the branches of 
continuous spectrum, see results in [HE]- This eigenvalue has bifurcated at p = from the lower 
end point ( A_ ) of the continuous spectrum. Our algorithm can follow this discrete eigenvalue 
on the imaginary axis, as illustrated in Figure ^p. The edges of the two continuous branches 
corresponding to positive and negative energy states in the massive Thirring model (|1.2|) are 



located on the imaginary axis at A+ = + 1) (upper branch) and at X_ = i{l — (lower 
branch) (see Section [2.3(1 . Here the discrete eigenvalue which lies on the imaginary axis is clearly 
below the edge of the lower branch ( A_ ) of the continuous spectrum. In the gap between the 
branches of continuous spectrum, the Evans function is still analytic, since there is still a 2:2 split 
of the dimensions of the stable and unstable manifolds. 




Figure 3: The real versus imaginary parts of the Evans function E{X) for the non-integrable case 
p = 0.1, V = 0.9 and 9 = 0.47r. The spectral parameter A varies parallel to the imaginary axis 
with a small offset off = 10~^^ . Going round clockwise from the upper left corner: 
(a): Overview of the Evans function. 

(h): First zoom into the area near E = 0. The colours are used to help identifying the lines in 
the next picture. 

(c) : Next zoom into this area. This picture shows that the two little loops cross right of the zero 
point. For visualisation purposes, we used off = 5 • 10~^ . off = the crossing of the loops is 
located exactly at the origin, indicating the discrete eigenvalue on the imaginary axis. The slight 
offset off = 5 • 10^^ moves the crossing of the loops to the right of the origin. 

(d) : Final zoom into the area near E = (note how small the scale is). We have increased the 
offset to off = 5 • 10^'^ to avoid problems with the smallness of the Evans function near A = . 

(e) : A topologically equivalent sketch of the Evans function. 

For < < 7r/2, the solitary waves are known to be stable even for non-zero p [HE]- In 
Figure 131 we show the Evans function for such a stable case. To determine the winding number, 
we need to zoom into the neighbourhood of E = 0. The behaviour near E = (Figure |3b~ 
d) and the global behaviour of the Evans function (Figure |3K) convey that the Evans function 
is topologically equivalent to a loop which does not contain the origin (Figure Efc). Hence the 
winding number is zero, confirming stability. 

For Tr/2 < 6 < IT , instabilities may arise for non-zero p. These instabilities are of an oscillatory 
nature in the sense that discrete eigenvalues leave the imaginary axis into the complex plane [Hll5j. 



Figure 4: The onset of the primary instability for p = 0.1 and V = 0.9 (the offset off = Oj. 
Plotted are the real part (continuous red line) and imaginary part (dashed blue line) of E{\) 
versus A. The vertical dotted line indicates the edge of the lower branch (\-) of the continuous 
spectrum. In the left picture with 9 = 0.502687r, the eigenvalue is still located on the imaginary 
axis. In the middle picture with 6 = 0.502707r, the eigenvalue has just merged with the end point 
of the lower branch (\-) of the continuous spectrum, giving a zero of the Evans function on the 
imaginary axis at the end of the continuous spectrum. In the right picture with 6 = 0.502727r, 
there is no eigenvalue on the imaginary axis anymore. 

The onset of the first instabihty is when the discrete eigenvalue in the gap of the continuous 
spectrum merges with the edge of the lower branch ( A_ ) of the continuous spectrum. In the 
6-p bifurcation plane, the instability curve starts at p = and 9 = 0.5 initially proportional to 
[H I15j. In Figure 131 the onset of the instability is illustrated at p = 0.1 and V = 0.9. At 
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Figure 5: (a): The real versus imaginary parts of the Evans function E{X) for p = 0.1, V = 0.9 
and 9 = O.Gvr. The spectral parameter A varies parallel to the imaginary axis with off = 5 - 10""^^ . 
(h): The Evans function in Figure\^ is topologically equivalent to this sketch of E{X). The 
winding number is clearly 2 , confirming that there is one pair of unstable eigenvalues for these 
parameter values. 

the bifurcation point the eigenvalue detaches from the imaginary axis at the edge of the lower 
branch ( A_ ) of the continuous spectrum and leaves into the complex plane. This can be explored 
by looking at the winding number of the Evans function. Figure El shows the real and imaginary 
parts of the Evans function when the instability has well occurred (at p = 0.1, V = 0.9 and 
9 = 0.67r). The winding number is 2, confirming that a pair of unstable eigenvalues is present. 

If we now decrease p (with 9 fixed) from the above described primary bifurcation of oscillatory 
instability, the primary unstable eigenvalue will merge with the edge of the upper branch ( A+ ) of 
the continuous spectrum. This is illustrated in Figure El If p is decreased below this bifurcation, 
the solitary wave is stable again. We see clearly the collision of the eigenvalue with the edge 



of the upper branch (A+) of the continuous spectrum. In 15. , a shghtly different perturbation 
is studied and it is proved analytically that in this case two bifurcation curves originate from 
(6, p) = (O.Svr, 0) . However, it is not clear what the relation with the branches of continuous 
spectrum is away from the bifurcation point. Our calculations suggest that for our perturbation 
there are again two bifurcation curves, one curve is related to a bifurcation from the edge of 
the lower branch (A_) and the other curve to a bifurcation from the edge of the upper branch 
(A-|_) and the primary unstable eigenvalue goes between these two branches. This seems a new 
observation, which has not been recognised earlier. 




Figure 6: Real (continuous red line) and imaginary part (dashed blue line) of the Evans function 
E{X) as a function of the spectral parameter A for V = 0.9 and 9 = O.Gvr (the offset off = 
5-W~^'^ ). (a): Above the bifurcation point at p = 0.002. (b): At the bifurcation point p = 0.0001 . 

Figures 0] and ini illustrate that our algorithm and the Evans function allow us to exactly follow 
the eigenvalues in the course of the bifurcation. By looking at the real and imaginary part of the 
Evans function the value of lambda which corresponds to a root of the Evans function can be 
located. This allowed us to answer the question from where on the imaginary axis the eigenvalues 
detach for the primary instability. 

It has been found analytically in Q that secondary bifurcations can occur for suitable values 
of p and 9 (we assume V fixed here for simplicity). Here a second pair of eigenvalues detach 
from the continuous spectrum into the complex plane. Again, analyticity of the Evans function 
allows us to compute the number of unstable eigenvalues by determing the winding number. In 
Figure 13 we show that the winding number equals 4 when p = 0.1, V = 0.9 and 9 = O.dn , 
indicating that a new pair of eigenvalues have detached. The behaviour near A = close to the 
origin of the Evans function behaves is similar to the one depicted in Figure |211. The onset of 
this secondary instability at p = OA and V = 0.9 is shown in Figure |H1 Figure |H1 shows the 
behaviour of the Evans function shortly before, at, and shortly after the secondary instability. 
In order to see the secondary instability and the emergence of the eigenvalue, one needs to look 
at the behaviour near E = 0. Note that the secondary instability is hard to observe in the full 
Evans function (left pictures in Figure |H1). In the middle picture, the Evans function has a zero, 
just after the edge of the upper branch (A+) of the continuous spectrum, but not at the edge of 
this branch. Hence for p > 0, the eigenvalue comes really from within the continuous spectrum, 
not out of one the endpoints of the branches of the continuous spectrum. Similar behaviour has 
been described by Sandstede & Scheel [23j . 

The results in this section show that our algorithm using the Evans function approach exer- 
cised on the wedge space provides a good and reliable diagnostic tool for the accurate detection 



Figure 7: (a): The real versus imaginary parts of the Evans function E{X) for p = 0.1, V = 0.9 
and 9 = O.Qvr. The spectral parameter A varies parallel to the imaginary axis with off = 5 - 10"-^^ . 
(b): The Evans function in (a) is topologically equivalent to this picture. The winding number is 
clearly 4, confirming that there are two pairs of unstable eigenvalues for these parameter values. 

of oscillatory instabilities. The analyticity of the Evans function allows for detection of complex 
eigenvalues by assuring the necessary conditions for the application of Cauchy's principle. Our 
algorithm is able to follow the location of the discrete eigenvalues in the course of their bifur- 
cations. The Evans function does not involve a discretization of the spectrum of the eigenvalue 
problem and hence does not fracture the continuous spectrum. It is therefore free of spurious 
unstable eigenvalues. Moreover, the exterior algebra is a platform for numerically stable shooting. 
As already stated, the method is not restricted to the particular coupled mode model H1.2|) but 
can be applied to other problems as well. 
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Figure 8: The real versus imaginary parts of the Evans function E{X) for p = 0.1, V = 0.9 
and 6 varying. The spectral parameter A varies on the imaginary axis (oS = 0). The green line 
helps to locate the emergence of the second eigenvalue. 

(a) : Just before the onset of the secondary instability at 9 = 0.790. The right picture zooms in 
at the neighbourhood of E = 0. 

(b) : At the onset of instability at 9 = 0.791. The right picture shows the real (continuous red 
line) and imaginary (dashed blue line) parts of the Evans function as function of A. This shows 
clearly that the zero of the Evans function occurs just after the edge of the upper branch (\+) of 
continuous spectrum (i.e., the cusp), hence inside the continuous spectrum and not at the edge. 

(c) : Just after the onset of the secondary instability at 9 = 0.792. The right picture zooms into 
the neighbourhood of E = . The loops with the green line now include the origin (E = 0), 
increasing the winding number by 2 and illustrating that an additional instability has occurred. 



